Comparison between two simple numerical models for the magnetoelectric interaction 

in multiferroics 
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We developed numerical calculations to simulate the magnetoelectric coupling in multiferroic 
compounds, using the Monte Carlo technique. Two simple models were used to simulate the com- 
pounds. In the first one, the magnetic ions are represented by a spin 1/2 2D Ising lattice of ions, 
and the electric lattice by classical moments, coupled one to one with the magnetic moments. The 
coupling between both lattices allows to the leading lattice, that is, the magnetic one, to change 
the orientation of the electrical dipoles in one direction perpendicular to the magnetic dipoles. This 
direction was chosen to accomplish the symmetry requirements of the magnetoelectric effect. In 
the second case, the magnetic lattice is also a 2D Ising lattice, but the electric momenta are in a 
lattice that also behaves as an Ising lattice, perpendicular to the magnetic moments. In this case, 
the one-to-one coupling of the electric and magnetic momenta is represented by a two- valued energy 
parameter, allowing the possibility of independent transition temperatures for both lattices. Both 
models contain three independent parameters. We studied the physical properties obtained with 
both models, as functions of the ratio of the three parameters. The results in both cases allowed 
us to compare changes in the physics of the models, and with the physics of compounds measured 
experimentally. 

PACS numbers: 75.10.-b, 75.10.Hk 
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Recently, there has been a revived interest in the re- 
search of multiferroics due to several new discoveries, and 
the possibility to use them technologically."'^ The electric 
and magnetic transitions are not necessarily correlated, 
but when it occurs - and the so called magnetoelectric ef- 
fect appears - the materials suggest possible use as mem- 
ories, etc. Besides the technical applications, several fam- 
ilies of those compounds present very rich physics. Just 
as an example, one of this compounds, LiNiP04, shows 
a phase transition where the electric lattice is not only 
first order, but in a very short range of temperature sev- 
eral incommensurable transitions appear.^ Even in front 
of those interesting phenomena, we found very few the- 
oretical papers on this subject, and most of them using 
very elaborated theoretical methods.'^ 

With the motivation presented above, we developed 
two simple numerical models, based on the Monte Carlo 
method, and the Metropolis minimization of total energy. 

We present here two simple and understandable models 
for the magnetoelectricity. The electromagnetic Hamil- 
tonian is solved for very simple cases, and the solutions 
present similitude with the reported experiments. The 
first model study the phase transitions at the same tem- 
perature, independently of the temperature value. The 
second model, a little more elaborated, allows the fer- 
roelectric and ferromagnetic transitions to occur at dif- 
ferent temperatures, and we studied the behavior of the 
model as function of this temperature difference. We 
compare the results between the models and with exper- 
imental cases. 



The physics behind the magnetoelectric effect consists, 
as seen in a bird's eye view, in the creation or orientation 
of electric dipoles by the magnetic moments or vice- versa. 
In the first case, the magnetic dipoles, which are perma- 
nent, when change their orientation, modify the lattice 
in such a way that the negative electric charges displace 
relatively to the positive. This is accomplished via spin- 
orbit coupling of the spins, changing the total energy 
as the orbit lattice Hamiltonian explains. Simplifying 
the model, the spin-orbit-lattice energy is calculated as 
a spin-lattice Hamiltonian. We simplified the calculation 
even more, representing the magnetic lattice as a 2D Ising 
lattice in all the cases. This corresponds with many real 
compounds, as the olivines mentioned above, where the 
structure of the real lattice presents separated planes of 
magnetic ions.^~® 

We assume that the magnetoelectric system is a set of 
magnetic dipoles, coupled via the exchange interaction, 
in a lattice with a distribution of electric charges, suscep- 
tible to change when the magnetic dipoles change their 
orientations. The change in orientation of the magnetic 
dipoles modify their environment, via spin-orbit interac- 
tion, creating local strains, and creating or orienting a 
set of electric dipoles in the lattice. We assume that 
our crystal suffers the strain in such a way that electric 
dipoles are oriented to a particular direction when the 
magnetic dipoles relax. 

The model Hamiltonian used in the models follows: 

H = Hm + He + Hme (1) 
where is the magnetic energy. He the electric energy 



2 



and Hme the magnetoelectric coupling. 

A. The first model 



The first approach for a solution of cq.(l) is obtained 
replacing the first term in the sum of the right side by 
a square sublattice of Ising magnetic moments, and the 
electric moments in the second and third terms by ran- 
dom oriented classical electrical dipoles, located in a sep- 
arated square sublattice. The Ising spins are coupled to 
their nearest neighbors only, and with periodic bound- 
ary conditions. The interaction Hamiltonian allows only 
nearest neighbors magnetoelectric interaction. Symme- 
try requires that the magnetic point group of the mag- 
netic moment is one of the 58 Shubnikov groups that 
allow magnctoclcctricity.^ This forces own magnetic mo- 
ments to have only one electric dipole as a nearest neigh- 
bor. The electromagnetic coupling is divided in two 
parts: the local interaction between the spin and the elec- 
tric dipole, and the lattice total electromagnetic energy, 
that takes into account the interaction between the elec- 
tric dipoles and their electric neighbors. As most of the 
electric parameters are measured perpendicularly to the 
magnetization,^'* wc chose the z direction for the mag- 
netic moments, and the x axis for the electric dipoles. 

The numerical solution of the problem was done using 
the importance sample Monte Carlo method, looking for 
the minimum in energy for our system. 

Thus, 

H =-J ^ (ji(7j - h^cji - /3 ^ PixPjx+l^Pix 

<i,j> i {ij} i 

(2) 

is the approximated Hamiltonian, where J is the ex- 
change coupling of the Ising magnetic spins ct, and P 
the electric dipoles. The symbol < i,j > indicates sums 
over the nearest neighbors only. The first and second 
terms constitute the magnetic energy, where we included 
the possibility of an applied or external magnetic field h. 
The third term represent the electric energy, proportional 
to the orientation of the neighboring electric momenta. 
As the system is to be ferroelectric, and the direction of 
the polarized dipoles the x axis of the crystal, we con- 
sidered only the energy coupling in that direction, which 
is represented by {i,j}, indicating sum over the two first 
neighbors located in the x axis. 

The interaction term, which represents a spin - lat- 
tice Hamiltonian, was separated into two parts. One of 
them is the local interaction between the spin and the x 
projection of the electric dipole, which makes that every 
transition of the spin changes simultaneously the dipole; 
the second, represented the last term in eq. (2) is the 
contribution of the lattice as a whole to the total en- 
ergy. As the local interaction is the same for every pair 



spin-dipole, we did not include it in the Hamiltonian. 
However, the meaning of this local part of the energy is 
important, as we will see below. 

1. Ferromagnetic case. 



Our first calculation was performed in a 100x100 2D 
lattice of Ising ferromagnetic spins coupled to 100x100 
electric dipoles, located in another square lattice, paral- 
lel to the magnetic one, and slightly shifted from it. The 
electric dipoles were oriented at random, together with 
magnetic lattice, to begin with infinite temperature. The 
temperature was then fixed to a value, and a Monte Carlo 
program, where the transitions are allowed following the 
Metropolis technique,^ is iterated the time necessary to 
obtain thermal eqiiilibrium of the system. Then, the re- 
sults arc used as the initial condition for the following 
temperature. The calculation was performed reducing 
the temperature in each step. 

The complete calculation was performed after a study 
of convergence in our model. As first step, we unconsid- 
ered the electric and interaction energies when we looked 
for the minimum. This means that the model is just a 
2D Ising system, moving electric dipoles together, and 
as expected, the magnetization follows the Ising model. 
The exact calculation published by Onsager allows a very 
good comparison, and the electric dipoles also feel a tran- 
sition at the same temperature. This calculation decided 
the size of the set of moments, which we selected as 
100x100 on this basis. 

The following step was to study the convergence when 
the parameters f3 and 7 in the model are different from 
zero. It is well known that the Ising model converges 
slowly near by the transition temperature, due to fiuc- 
tuations, and the equal value for the energy when the 
system is oriented in any of the both possible directions. 
This is easily solved with the addition of the small ex- 
ternal field h; however, we observed that for particiilar 
values of 7, the convergency is the slowest. This can be 
explained by the fact that 7 appears in the Hamiltonian 
as an extra external field, in some manner. We can use 
the local coupling of the spin-dipole pair to write 

Pix = Pix\o-i\ = ai\Pix\ (3) 

for the modulus of cTj is always the unity; that can be used 
to write the second and last terms in the Hamiltonian as 



l^Pix-hY,(^i = -^{h--f\Pix\)<Ti (4) 

i 

which shows that the value of 7|Pia;| appears as an ex- 
tra magnetic field. The mean value of \Pix\ annulate 
the external field when 7/ J « 0.02 for an applied field 
h/J = 0.01, and the convergency is the slowest for this 
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value of 7. Taking this into account, we found that it 
was necessary 5000 iterations per spin to obtain thermal 
equilibrium and other 5000 to get the mean values of 
the energy and magnetization; this numbers were used 
in every case. 

Eq. (4) may also be used to analyze the meaning of 
the /? parameter. The first and third terms in eq.(2) can 
be written 



= -J ^ gjo-j - /3 \PixWi\P]xWj 

<i.j> {i,j} 
= -^iJ + PiPixWPjxDcriCrj ~ J ^ O-jCTj 

So it can be seen that the /3 parameter makes the ex- 
change anisotropic, modifiyng its value in the x direction, 
leaving the coupling in the y direction unaltered. 

We performed the calculation as function of the tem- 
perature for different values of the /? parameter for 7 = 0; 
then repeated the calculation for /3 = to study the de- 
pendence of the results from 7, and finally we made a 
complete study of the form and transition temperature 
of the system as a function of both parameters. As J de- 
termines the transition temperature for the Ising model, 
we used it as unit of energy for the whole system. 

To make the results clear, we made calculations as 
function of the temperature for different values of /3, 
when 7 is zero - meaning that the electric interaction 
is bigger than the magnetoelectric one. After that, we 
calculated the minimum as function of 7, when /3 is zero. 
The complete calculation, with both parameters different 
from zero gives a clear vision of the total behavior of our 
model. The results for the ferromagnetic case are shown 
in Figs. 1, 2, and 3. The first and second figures show 
the change in the shape of the transition caused by the 
/? and 7 parameters. The effect of the /3 parameter is 
principally to shift the transition, as seen in Fig. 2. The 
7 value can be positive and negative, and the effect is 
to broaden the transition, and invert the magnetization- 
polarization when the sign of it is changed. Fig. 3 shows 
the complete dependence of the transition temperature 
depending of both parameters. 

Figs. 1, 2 and 3 show the results for the ferromagnetic 
case. Fig. 1 shows the effect of the 7 parameter: as the 
transition is determined by J, the effect of 7 is to broaden 
the transition. Fig. 2 shows the effect of f3, which is to 
shift the transition without great modifications in the 
shape of it. Fig. 3 resumes the complete model, show- 
ing the effect of both parameters together. The induced 
electric polarization appears at the same temperature as 
the magnetic transition in all the cases, as it is defined 
by the model. As the model only allows, both P/Pq and 
M/Mq curves coincide. 
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FIG. 1: (color online) Normalized magnetization (and electric 
polarization) as function of the temperature when P — in 
the first model. 
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FIG. 2: (color onfine) The same temperature dependencies as 
in Fig. 1, when 7 = 0, in the first model. 



2. AntiferromagneUc case 

The antiferromagnetic case was treated similarly. It 
requires a negative value of J, but several changes in 
the other terms of the Hamiltonian are necessary. The 
magnetization of both sublattices will be coupled to the 
electric lattice opposed, in order to obtain the required 
ferroelectricity. The magnetic field is set to zero, because 
only can be directed parallel to one of the magnetic sub- 
lattices. The convergence is slowest for zero field, which 
we used to obtain the values for convergence. Our re- 
sults in this case are very similar to those above, and 
we lack here of space to show them completely. Fig. 4 
presents the transition temperature dependence for this 
case, which can be compared with the ferromagnet. The 
complete results will be publish elsewhere, together with 
a more sophisticated model for the spin-lattice coupling. 
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FIG. 3: (color online) The general results for the first model. 
The transition temperature is the same for both the electric 
and the magnetic transitions (see text). 
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FIG. 4: (color online) The general results for the first model 
in the antiferromagnetic case. The transition temperature is 
the same for both the electric and the magnetic transitions 
(see text). 



FIG. 5: (color online) Normalized magnetization and elec- 
tric polarization as functions of T for /3/J = 1 in the second 
model. 



B. The second model 

As seen above, the first model does not contain the 
capacity to allow different temperatures for the electric 



and the magnetic transitions, giving us only the changes 
generated in the shape of the transitions by the magneto- 
electric coupling. We decided to develop a second model, 
where the transition temperatures are independent. To 
maintain the simplicity of the model, and the number 
of independent parameters reduced to three, we included 
an energy A — e2 — Si for the pair of magnetic-electric 
moments. When the spin is up, and the electric dipole 
points the left, or when the spin points down, and the 
electric dipole to the right, they will have an energy A 
higher than in the other two cases. If we make A ^ cxo 
we recover the first model again, for the system will be 
in the lower state all time. 

We changed some other things in the new model. In- 
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FIG. 6: (color online) Calculated magnetization and polariza- 
tion per spin/dipole, for the second model when P/ J — 2 as 
functions of T. It can be seen that M/Mo is strongly distorted, 
while P/Po is not. 



FIG. 7: (color online) Results for /3/J = 0.5 in the second 
model. In this case, the magnetic function is not distorted, 
but the polarization is. 



stead of the classical electric dipoles oriented at random 
at T — oo, we substituted the electric lattice by an Ising 
lattice, oriented in the x direction, that is, as formerly, 
the electric dipoles form the ferroelectric part of the lat- 
tice when they are oriented in the positive x direction. 
We excluded the local interaction between the pairs, con- 
sidering that the electric interaction is only in the x di- 
rection, let's say, the tails of the interaction after first 
neighbors cutted off. The result is that we substitute 
the magnetoelectric interaction for this two-level system. 



Hence, the complete Hamiltonian is, in this case 

H = -j o,Gj-hY,<y]-p ("^^ 

<i,j> i <i,j> i 

where the symbols are the same as before, and the Si the 
energy of the pair magnetic-electric momenta as defined 
above. We used again the exchange coupling parameter 
J as energy unit; as can be seen , this model has two in- 
dependent transition temperatures for the magnetic and 
electric lattices, as J and /3 are independent in this model. 

To use the Monte Carlo method again, we need to 
preserve the mathematical requirements for it, then our 
model needs to behave as a Markovian one. To accom- 
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FIG. 8: (color online) Results of the second model when 
A/ J = 0.5 and changing the /3 parameter for P/J < 1 . 
The shape of the polarization curves are distorted and shifts 
as the electric transition temperature is increased. 



plish this requirement, the minimum of energy is calcu- 
lated through the following steps: 

a - As in the first model, the initial state of the system 
is at T — ^ oo. Then the value of the temperature is 
inserted in the calculation. 

b - One number A is then chosen at random between 
zero and one. If this number is zero, we invert the cor- 
responding spin; instead, for A = 1, the inversion is per- 
formed on the electric dipole. 

c - One of the pair of moments is chosen at random. 
If A = 0, we calculate the energy difference if the spin is 



FIG. 9: (color online) Calculated polarization and magneti- 
zation for A/ J — 0.5, for values of the /? parameter when 
P/J ^ 1. 



inverted. From eq. (6), this energy difference will be 

AEi^2<j,{JSs + h) + <j,P„A (7) 

where ai is the chosen spin, Ss the sum of the first spins 
which are nearest neighbors to it, and Pix the component 
of the electric dipole of the pair. If AEi is negative, the 
spin is inverted. If AEi is positive, we use the Metropolis 
comparison: if a new random number < r < 1, when 
compared with the energy population is less than that 
(that is r < exp (— Ai^i/fesT)) the spin is inverted too. 
If not, the spin is left unaltered. 

d - If A = 1 the electric dipole is inverted; the energy 
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FIG. 10: (color online) Electric and magnetic transition tem- 
peratures, for A/J = 0.5 as depending of the /? parameter. 



1. Convergence. 

Differently from the first case, this second one does not 
contain a strong local coupling for the pair, and even the 
local change is not always decided by the spin system. It 
was necessary, then, to realize an independent study for 
every pair of the parameters, A and /?. We do not believe 
that the insertion of the whole convergence study could 
be interesting for the reader, so we just mention that the 
number of Monte Carlo steps per spin for convergency 
varies from 5 to 50 thousand steps, the small value for 
A/ J = 0, ^/ J = 1 and the maximum for A/ J = 1, /3/ J = 
2. 



2. Results 



1.0 



0.8 



0.6 



S 0.4 



0.2 



Second model 
First model 



10 



12 



kT/J 



FIG. 11: (color online) The magnetoelectric coefficient as cal- 
culated for our models. As expected, they tend to zero at 
T=0, as the magnetizations and polarizations saturate. 



difference in this case is 



Ai;2 = 2P,,PSd + (T^P^x^ 



(8) 



Here the value of is the sum over electric dipoles which 
are nearest neighbors to the one chosen. The other sym- 
bols have the meaning above. Again, if IS.E2 is negative, 
the dipole is inverted, and if positive, the random number 
r is calculated and used to compare with the populations 
in order to decide the inversion of the dipole. 

e - The procedure from b) to d) is repeated as many 
times as necessary to obtain convergence to thermal equi- 
librium as defined for the first model. Then the temper- 
ature is reduced further and the calculation repeated to 
equilibrium. 

The model respects the symmetry requirements for 
magnetoelectricity. The system does not contains tempo- 
ral nor spacial inversion, so this exigency is accomplished 



We performed complete calculations for several cases, 
described below 

1 - Thermal dependence of polarization and magneti- 
zation when j3 / J — 1 (that is, the transitions occur at 
the same temperature) and A varies. 

2- The same study, for /?/ J > 1. 

3 - The same again, when (3/ J < 1. 

4 - The thermal dependence, as function of /3, for 
A/J = 0.5. 

Now we will describe the results for every case from 1 
to 4. 

1 - As the electric dipoles are not strongly coupled to 
his magnetic neighbor, and they can assume both ori- 
entations, the shapes of the transitions are different, as 
can be seen in Fig. (5). As there is not an applied elec- 
tric field, the relative orientation of the polarization and 
the magnetization could be at random, as can be seen in 
the figure. The transition temperature increases with the 
value of A, meaning that the coupling helps to mantain 
the system ordered. 

2 - Fig. (6) shows the results for /3/J ~ 2. As can 
be seen, the transitions differ strongly in shape, and the 
curve of magnetization is shifted to accompany the po- 
larization transition, as A increases. This coincides with 
the idea that the second model limits with the first, when 
A — cx). Anyway, we were not able to calculate this case 
for high values of A, for the time of convergence increases 
too much. 

3 - The case where the electric transition occurs at 
lower temperature than the magnetic one is presented in 
Fig. (7). Here, the magnetic transition is not deformed, 
and the electric one is distorted, with the tendency to ac- 
company the magnetization for A — 00. Again, we were 
limited to values of A allowed by the time of convergency 
of the calculation. 

4 - The calculation was performed for A/J — 0.5 for 
different values of the transition temperatures ( J ^ /3). 
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The results are shown graphically in Figs. (8) and (9). 
We separated the results when the electric transition is 
at lower temperature than the magnetic (Fig. (8)), and 
the opposite, when /3 > J (Fig. (9)). It can be observed 
in both cases that the transition occuring at higher 
temperature remains almost undistorted, and the; one 
whose transition temperature is smaller, distorced and 
shifts. Fig. (10) shows the transition temperatures, elet- 
ric and magnetic, when A/J = 0.5 as functions of /?. Of 
course, both curves equal when (3/ J = 1. The magnetic 
transition tends to saturate for small or great values 
of j3, while the electric transitions behave almost linearly. 



III. CONCLUSIONS AND PERSPECTIVES. 

Magnetoelectricity and magnetoferroics are studied ex- 
perimentally using the phenomenological free energy ob- 
tained just from symmetry and the (possible) interaction 
between magnetic and electric fields, as follows^: 

F{t,t) = Fo-PfEi-MfHi- 

~2^o^ijEiEj — -/iQii^jHiHj — aijEiHj — ... 

and we have: 

dF „ 

Pi = -Q^,=Pi +^o^i3Ej+aijHj + ... 

dF a 
" -^=^i +MoM^i^^.• + ai^^^J•-••• 
where it can be seen that the experimental measurement 
of the magnetoelectric tensor, a^ , is performed looking 

for the difference of the observed magnetization (polar- 
ization) with and without an external magnetic (electric) 



field. To calculate that parameter, we followed the same 
procedure, calculating for every temperature the mag- 
netic polarization with and without an applied field, that 
is 

P - P^ 
a,,(T)^^^^ 

as the experiments are made.^ 

Fig. (11) presents the results for both models, when 
the transition temperatures are the same (J = /3) in the 
second model. The magnetoelectric coefficients go to 
zero at T = 0, which is expected in a system that sat- 
urates magnetically and electrically too. Our models do 
not include the possibility to change the energy saturated 
at K, but the experiments (see ref. 3) show a remanes- 
cent value of the parameter. Our model does not include 
any effect in other than the x axis, thus, we only obtained 
the axz magnetoelectric coefficient within both models. 

The transition in ref. 3 is first order, that is, the coeffi- 
cient docs not exist for temperatures above the magnetic 
transition, contrary to our models, where both transi- 
tions are second order, and as that, the curves in Fig. 
(11) extend to high temperatures. 

Resuming, our calculation arrives to many similitudes 
and differences with experiment. We believe that this 
can be the simpler way to simulate real systems, and 
developing more elaborated spin - lattice terms in the 
Hamiltonians will help to interprete the experimental re- 
sults in an easy way. 
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